Chapter 8 Captive-bred vs wild

load("data/beta_filtered_capwild.Rdata")

wild_pre_samples <- sample_metadata %>% 
  filter(diet!="Post_grass") %>%
  dplyr::select(sample) %>% pull()

genome_counts<- genome_counts %>% 
  column_to_rownames("genome") %>% 
  select(all_of(wild_pre_samples))%>% 
  rownames_to_column("genome")

sample_metadata <- sample_metadata %>% 
  filter(diet!="Post_grass")

beta_q0n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1)

genome_counts <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_tree <- keep.tip(genome_tree, tip=genome_counts$genome)

beta_q1p <- genome_counts %>%
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, tree = genome_tree)

genome_counts_filt <- genome_counts[genome_counts$genome %in% rownames(genome_gifts),]
genome_counts_filt <- genome_counts_filt %>%
  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0))%>%
  rownames_to_column(., "genome")

genome_gifts1 <- genome_gifts[rownames(genome_gifts) %in% genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]

dist <- genome_gifts1 %>%
  to.elements(., GIFT_db) %>%
  traits2dist(., method = "gower")

beta_q1f <- genome_counts_filt %>%
#  remove_rownames() %>% 
  column_to_rownames(., "genome") %>%
  filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
  select_if(~!all(. == 0)) %>%
  hillpair(., q = 1, dist = dist)

8.1 Permanova analysis

set.seed(1234)

Richness

sample_metadata_row<- column_to_rownames(sample_metadata, "sample") 
sample_metadata_row <- sample_metadata_row[labels(beta_q0n$S), ]

betadisper(beta_q0n$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.004224 0.0042240 7.8055    999  0.003 **
Residuals 22 0.011905 0.0005412                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.008
Wild       0.010583      

adonis2(beta_q0n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q0n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.01526961	0.2489907	7.29391	0.001
Residual	22	0.04605641	0.7510093	NA	NA
Total	23	0.06132602	1.0000000	NA	NA

Neutral

betadisper(beta_q1n$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.006925 0.0069254 0.5404    999  0.491
Residuals 22 0.281957 0.0128162                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.476
Wild        0.47005      

adonis2(beta_q1n$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1n$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	1.538238	0.2539205	7.487472	0.001
Residual	22	4.519715	0.7460795	NA	NA
Total	23	6.057953	1.0000000	NA	NA

Phylogenetic

betadisper(beta_q1p$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.007671 0.0076708 1.8217    999  0.179
Residuals 22 0.092637 0.0042108                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.187
Wild        0.19084      

adonis2(beta_q1p$S ~ diet, 
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1p$S))), 
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.3720264	0.3082818	9.804858	0.001
Residual	22	0.8347476	0.6917182	NA	NA
Total	23	1.2067741	1.0000000	NA	NA

Functional

betadisper(beta_q1f$S, sample_metadata$diet) %>% permutest(., pairwise = TRUE) 

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.003070 0.0030698 1.7361    999  0.171
Residuals 22 0.038901 0.0017682                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Pre_grass  Wild
Pre_grass           0.176
Wild        0.20119      

adonis2(beta_q1f$S ~ diet,
        data = sample_metadata %>% arrange(match(sample,labels(beta_q1f$S))),
        permutations = 999) %>%
        broom::tidy() %>%
        tt()

term	df	SumOfSqs	R2	statistic	p.value
diet	1	0.01079711	0.187686	5.083124	0.103
Residual	22	0.04673041	0.812314	NA	NA
Total	23	0.05752752	1.000000	NA	NA

8.2 Beta diversity plots

8.2.1 Richness

beta_q0n$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Post_grass" = "Captive-bred grass diet", "Wild" = "Wild"))+
    #   stat_ellipse(aes(color = beta_q1n_nmds$Groups))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(face = "bold", size = 12),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 12),
      legend.title = element_text(size = 14),
      legend.position = "right", legend.box = "vertical"
    ) +labs(color='Origin')

8.2.2 Neutral

8.2.3 Phylogenetic

beta_q1p$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
    geom_point(size = 4) +
    scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
    geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
    theme_classic() +
    theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')

8.2.4 Functional

beta_q1f$S %>%
  vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  dplyr::left_join(sample_metadata, by = join_by(sample == sample)) %>%
  group_by(diet) %>%
  mutate(x_cen = median(NMDS1, na.rm = TRUE)) %>%
  mutate(y_cen = median(NMDS2, na.rm = TRUE)) %>%
  ungroup() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = diet)) +
  geom_point(size = 4) +
  scale_color_manual(values = diet_colors,labels=c("Pre_grass" = "Captive-bred no grass diet", "Wild" = "Wild"))+
  geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
  theme_classic() +
  theme(
      axis.text.x = element_text(size = 12),
      axis.text.y = element_text(size = 12),
      axis.title = element_text(size = 20, face = "bold"),
      axis.text = element_text(face = "bold", size = 18),
      panel.background = element_blank(),
      axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
      legend.text = element_text(size = 16),
      legend.title = element_text(size = 18),
      legend.position = "right", legend.box = "vertical"
    ) +
    labs(color='Origin')